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A kinetic one-dimensional Ising model is coupled to two heat baths, such that 
spins at even (odd) lattice sites experience a temperature T e (T ). Spin flips occur 
with Glauber-type rates generalised to the case of two temperatures. Driven by the 
temperature differential, the spin chain settles into a non-equilibrium steady state 
which corresponds to the stationary solution of a master equation. We construct a 
perturbation expansion of this master equation in terms of the temperature difference 
and compute explicitly the first two corrections to the equilibrium Boltzmann distribu- 
tion. The key result is the emergence of additional spin operators in the steady state, 
increasing in spatial range and order of spin products. We comment on the violation 
of detailed balance and entropy production in the steady state. 



I. INTRODUCTION 



In recent years, considerable interest has focused on non-equilibrium stationary states (NESS) in 
interacting many-body systems, induced by deterministic or stochastic microscopic dynamics QJ^]. 
In the stochastic case, the starting point is typically a master equation, i.e., a continuity equation 
for the (time-dependent) configurational probabilities. Different models are characterised by their 
microscopic transition rates. Both, equilibrium and non-equilibrium stationary states, have time- 
independent macroscopic observables. However, only in the equilibrium case it is possible to compute 
these observables without explicit reference to the imposed dynamics, in the framework of Gibbs 
ensembles. In stark contrast, NESS and their properties depend generically on the details of the 
dynamic rules. Since a unifying theoretical description of NESS is still lacking, most progress to date 
is made by studying specific models. 

In some of the simplest models, a NESS is established by driving the system with an external field. A 
characteristic example is the (fully periodic) Ising lattice gas, subject to a uniform electric field which 
induces a nonzero particle current A second class of models involves two temperature baths. These 
can be coupled either to the system boundaries or ac t throughout the bulk: in the latter case, 

each temperature bath controls a (translationally invariant) subset of configurational transitions. For 
instance, in an Ising lattice gas, particle-hole exchanges along certain spatial axes might be coupled 
to a higher temperature [HQ]. Or, in the non-conserved case, spin-flips on a selected sub-lattice might 
occur at a different temperature . In all of these cases, equilibrium can be recovered upon letting 
a key parameter vanish - e.g., the strength of the bias, or the temperature difference. In this spirit, 
these models allow us to probe the effect of non-equilibrium perturbations on equilibrium properties, 
and we may ask to what extent universal (i.e., long- wavelength, long-time) features are affected. 
One finds, generically, that non-equilibrium perturbations are far more relevant, in a renormalisation 
group sense, in systems with conservation laws 0Jl^,|ll|. In contrast, it can be shown that Ising-like 
models with non- conserved dynamics remain in the universality class of the equilibrium Ising model, 
even if the usual Zj£ -symmetry of the Ising model is broken \L^ . While these results are clearly 
of major theoretical interest, they do not reveal how NESS configurational probabilities differ from 
their equilibrium counterparts: How, for instance, is the Boltzmann distribution for the equilibrium 
Ising model modified, when a second, different temperature bath is coupled to some spin-flips or - 
exchanges? Clearly, such modifications must be severe if universal behaviour is to be affected. Yet, 



'Present address: MPI fur Physik komplexer Systeme, Nothnitzer Str. 38, 01187 Dresden, Germany, e-mail: 
frank@mpipks-dresden.mpg.de 
^e-mail: schmittm@vt.edu 



even if long-wavelength, long-time properties remain effectively Ising-like, fundamental changes in the 
configurational probabilities should arise from the breaking of detailed balance. 

It is with these motivations, questions and expectations in mind, that we turn to the arguably 
simplest non-equilibrium Ising-type model, namely, a one-dimensional interacting spin chain with 
spin-flip dynamics. The rates are a simple but nontrivial generalisation of the well-studied Glauber 
p3[ | rates: spins at odd (even) sites are coupled to a temperature T (T e ). In dimensions d > 2, 
this model exhibits an order-disorder phase transition which should be in the Ising universality 
class, according to renormalisation group arguments fl2]| . Detailed Monte Carlo simulations in d = 2 
confirm this expectation. In d = 1, an exact solution for the two spin correlation functions has 
been found jli]]. In the following, we seek an analytic expression for the steady state probability 
distribution, i.e., the solution of the master equation, for this model. While we have not been able 
to find an exact solution, our perturbative analysis provides some insight into how the equilibrium 
distribution is modified when a (small) temperature difference is present. 

As energy is fed into a non-equilibrium stationary state from one reservoir, it must also be dis- 
sipated: in general, a NESS constantly produces "entropy" . Possible definitions of "entropy" and 
"entropy production" in non-equilibrium systems have been discussed intensely in recent years, from 

Tq . However, different perspectives 
have been developed rather independently, so that the relations between them still remain somewhat 
unclear. Therefore, we find it worthwhile to compare results from different approaches to entropy 
production, using our model. 

We conclude with two comments: First, while some analytic results for steady-state distributions 
are available, they are confined to two classes of systems: first, d = 1 lattice gas models, restricted to 
excluded volume interactions, such as the asymmetric exclusion process and its relatives |l7|-p^|, and 
second, very special (d = 1) spin systems whose master equations are solved by the Ising Boltzmann 
factor Second, even though our main emphasis rests on fundamental aspects, we note briefly 

that two-temperature systems of this type can be established in real systems, and maintained for 
significant amounts of time. For example, via nuclear magnetic resonance in an external magnetic 
field, a lattice of nuclei in a solid can be prepared to have a certain spin temperature pl[ |. Thus, in a 
crystal with two sub-lattices of different nuclei the temperatures of these sub-lattices can differ. 

This article is organised as follows. We begin by introducing our model and its master equation. The 
following section describes our perturbative approach towards finding the stationary distribution. We 
then compute the two leading orders, aiming to highlight some generic properties of the perturbation 
series. We next turn to the discussion of two more general questions, namely, first, how spin flips violate 
the detailed balance condition, and second, whether different definitions for the entropy production 
rate lead to the same results, for our model. We conclude with an outlook and some open questions. 
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II. THE MODEL 



Our model is defined on a one-dimensional spin chain, of length N (with N even). Each site i 
carries a spin variable, cr^, which can take the values ±1. We impose periodic boundary conditions, so 
that (7/v+i = Spins at even (odd) lattice sites are coupled to a heat bath at temperature T e (T a ). 
Thus, a configuration {er} = {cri, 02, ■ ■ ■ , cjv} evolves into a new configuration by flipping a randomly 
selected spin, e.g., spin <7j, with a rate ]14[ 

Wi(<Ti -> -<7i) = 1 - —ai(ai-i + a i+ i) , (1) 

where 

_ f 7e = tauh(2J/k B T e ) , i even , . 

7i ~ \ 7 = tanh(2 J/k B T ) , i odd . [ ' 

Here, J denotes a nearest-neighbour exchange coupling between spins. From now on, we use di- 
mensionless units for inverse temperature, i.e., (3 e = J/(/csT e ), etc. The rates respect the usual 
Z>!-symmetry of the Ising model. The full stochastic dynamics of the system can be cast in terms of 
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a master equation for the time-dependent configurational probability p({a};t) 



N 

d tP {{a}-t) = -► -<Ti)p{{<r};t)+Wi(-<n -> <n)p ({<j [i] }; i)] (3) 

Here, {crM} differs from {er} by a flip of the i-th spin. A trivial time scale has been set to unity. 

An alternate description of our spin chain is the usual defect (domain wall) picture. A unique 
defect configuration {a} is associated with a pair of spin configurations, namely {a} and its image 
under Zj*, {— a}. The defects are located on the dual lattice (i.e., the bonds), with a% = (1) if 
<7j(7j+i = 1 (— 1). A spin flip er, — * — Oi in the spin chain {<t} can give rise to three different types of 
processes in the defect state {a} depending on the three neighbouring spins {(Tj_i, Cj, <j i+ i}: (a) if 
<Ji_i = — Ui = <Ji+i, two adjacent defects annihilate each other, with rate 1 + j e if i is even (1 + 7 D if 
i is odd); (b) if <7,;_i = <Ji — (Ti+i, two adjacent defects are generated simultaneously, with rate 1 — j e 
if i is even (1 — 7 C if i is odd); and (c) if <7j_i = — Cj+i, a defect diffuses to a neighbouring lattice site, 
with rate 1. 

Clearly, for uniform temperature T = T e = T , this model reduces to the well-known Glauber 
dynamics |0|, with homogeneous rates w(ai — > — cr,). The associated steady-state solution is just the 
Boltzmann distribution, 

lim p({a}; t) = \ exp (-H/k B T) = \q Q ({a}) (4) 

t—>oo Z Zj 

where H is the Ising Hamiltonian for a spin chain, 

N 

H({<j}) = -JY / ^+i (5) 

i=l 

and Z is the (canonical) partition function. Clearly, H ({a}) just counts the number of defects, 
n({cr}) = X)i=i a ii m configuration {a} via J -1 ff ({cr}) = 2n({a}) — N. The equilibrium distribution 
<7o ({f}) /-Z will serve as the unperturbed reference solution for our perturbative analysis, to be outlined 
below. For later reference, we note that detailed balance holds in the equilibrium case, i.e., 

w(oi -> -<7i) q ({a}) = w(-ai -> <Ti) q ({cr W }^ (6) 

for any spin- flip. 

Returning to the case of two temperatures, it is easy to see that the (time continuous) Markov 
process, characterised by Eqn. (||) is ergodic, since every configuration {a} can be reached after a 
finite time from every other configuration {a 1 }. Hence, any initial condition will converge, for t — > oo, 
to the unique stationary solution, q({a}). The determination of q({cr}) is our goal here. Being 
the full configurational probability distribution, it contains all information about the system in its 
non-equilibrium steady state. 

Since our NESS violates detailed balance (cf. Section 4), the determination of q({a}) is highly 
non-trivial. As a starting point, we restrict ourselves to the regime where the two temperatures differ 
only slightly from each other, so that perturbative methods can be invoked. In particular, we are 
interested which new spin operators are induced by the coupling to two heat baths. Writing the 
stationary distribution as 

q = I exp (V ({*})) , (7) 

we seek the "potential" function V ({c}) in perturbation theory (Z is determined through the nor- 
malisation condition). 



III. PERTURBATION THEORY 

In order to set up a perturbative treatment of the master equation (S) , it is convenient to decompose 
the flip rates (□) into an equilibrium-like contribution and a non-equilibrium perturbation 
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Wi{(Ji 



1 -. <Ji{?i-l + CTi+lj 



lo - le 



(-iy 



Ui (<7 i — l + <Ji+l) 



w(<Ji-i, (Ji, (Ji+i) + ° e Ai(ai-i, <Ji, a t+ i). 



(8) 



The term in the first [...] bracket, w), is just the homogeneous Glauber rate, at an effective temper- 
ature defined via tanh(2/3) = 7 = (j e + j ) /2. The second term is explicitly proportional to the 
"temperature" difference, 



lo - le 



0) 



and hence captures the non-equilibrium aspect of our dynamics. Thus, d serves as a suitable expansion 
parameter for a perturbation theory near equilibrium. We note that it is restricted to the interval 
< \d\ < 0.5 where |d| = 0.5 corresponds to one temperature being infinite and the other zero. 
We seek the stationary solution q ({a}) of Eqn. (||), namely 



= ^2 [- M°j-ii + d A j( CT i-i! a v <i {{<?}) 

+ -<Tj, (Tj+l) + d Aj(aj-i, <7j, cr i+ l)) 9 ({o^ 1 } 

We assert [||| that g ({c}) can be written as a perturbation series, 

q({*}) = Z- 1 qo ({a}) (l + fyQ„(M)] , 



(10) 



(11) 



where go ({ <J }) = exp(/3 ^^tJifJi+i) is the Ising potential at inverse temperature /3. Inserting this 
ansatz into Eqn. (|l0), and ordering terms in powers of d, uniqueness demands that the coefficient 
accompanying each power d n , n = 1, 2, ... vanishes. The lowest order is satisfied by construction. 
Higher orders are now evaluated recursively: the coefficient of d n , for n > 1, can be written as the 
sum of two terms which must cancel: = Si + S^ 1 . Here, s[ n ^ contains contributions of the form 
AjQ„_i and is therefore explicitly calculable with the help of results from the (n — l)-th order. The 
second term is generated by applying the equilibrium rates w(. . . ) on the unknown Q n . Even though 
we cannot invert the Glauber Liouvillean w exactly, we can infer the structure of Q n , based on the 
knowledge of Si and the simplicity of w (which contains only nearest-neighbour spin operators). 
Of course, at higher orders, the structure of Si becomes increasingly complex, so that explicit 
computations become quite cumbersome. However, at lower orders, this program is quite feasible. In 
the following, we illustrate the process for the first-order correction, and provide a few pointers for 
the second order. 

At order d, we need to find Qi {{a}) such that 

N 

= [-^(^-I'^'^+i)*) (M) + Ajiaj-^-a^aj+^qo ({ct [j1 } 

3=1 
N 



3=1 
_ c(l) 



X! \- w { a 3-i> a j> a 3+i)<lo (W})Qi (M) + w((jj-i,-(jj,a j+ i)qo ({<T b] }j Qi ({cr [j1 } 

3 = 



s 



(1) 



(12) 



where Si and S^' represent the first and second sum over j, respectively. We first compute 

c (i) /r „ sinh(4/5) A 
#1 = -9o (M) 2 1^ (-1)^3^+2 , 

3 = 1 



,(D 



(1). 



(13) 
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which suggests the following ansatz for Qi({cr}) 



N 



h(M)=A^(-l)V i (7 i+2 , (14) 



where the (real) parameter A needs to be determined. Inserting (|l4j) into Sg , we obtain 

N 

Si 1] =4X q (-1)'' °3°3+* > 



whence 



A = - Sinh(4 ^. (15) 



Thus, we arrive at the first order correction to the stationary distribution: 

9 (M) = | 90 ( W) ^1 - d t^L^M. J2(-iy a> + 0(d 2 ) . (16) 

Hence, in first order perturbation theory the temperature difference induces an interaction between 
next-nearest neighbours with a sign difference for even/odd pairs. The correction has a simple, intuitive 
interpretation in the defect picture. If n e ({o~}) and n ({a}) denote the number of defect pairs centred 
on even or odd sites, respectively, it is easy to see that 



1 N 



i=l 

i. e., Qi ({<r}) is proportional to the difference between the number of "even" and "odd" defect 
pairs. For T e > T Q , d is positive, so that, e.g., a spin configuration with a single even defect pair 
acquires a higher statistical weight than a configuration with a single odd pair. In contrast, those 
two configurations are degenerate in equilibrium, i.e, they have the same weight. Clearly, the high 
degeneracy of the equilibrium probabilities is reduced by the perturbative correction: Now, not only 
the total number of defects matters, but also how defect pairs are distributed over the two sub-lattices. 

Turning to the second order contribution, the calculations, while becoming more involved, proceed 
in essentially the same manner. Collecting all terms of Eq. ( |l0| ) with a factor d 2 the condition for the 
correction Q 2 ({c}) reads 

JY 



= °J5 9o (M) Qi (M) + Aj(oj-i, -a-j, (Tj+i) 

3=1 

N 

x 90 ({a®}) Qi (V [J] }) ] + £ [-^3-1, o-i+i)) 90 (M) Q2 ({*}) 

3=1 

+ w(a^ u -a 3 , a j+1 ) q Q ({<!«}) Q 2 ({a®})] =: S (2) + S (2) , (17) 
(2) (2) 

where again the two terms S{ and S 2 are defined as the first and the second sum over j, respectively. 
We obtain 

S{ 2) =4X 2 q ({a}) <M ^(-1)^*^+2 j 

N - 1 

+4 Y + o-iCi+i a i+2 a i+3 ) - coth(4/3) (aia i+1 + aia i+3 )] > (18) 

»=i J 
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with A given in Eqn. @. A gain, a suitable ansatz, based on the structure of S[ , allows us to 
determine the second-order correction Q2 ({c}) 




N ~ \ 

-4 ^ [{aiPi+2 + cri(Ti + icri + 2cri + 3) - coth(2/3) {u. L a i+ i + <Ti<T l+3 )] > 

4=1 J 



(19) 



Remarkably, we now encounter both, higher- order spin operators and next-next-nearest neighbour 
interactions. 

Clearly, (PQ2 ({c}) should be added to the first order solution, Eqn (|l6[). It is instructive to write 
the stationary distribution in exponential form, via q ({er}) = Z^ 1 exp E^Lo to illustrate 

the differences from the equilibrium form even more succinctly. To 0(d ), we obtain: 

V Q ({a}) = ^a i a i+1 , ^i(W) = A^(-l)V if 7 i+2 , (20) 

i i 

V 2 ({cr}) = 2A 2 ^ [o-jcr i+2 + o \o - i+1 a i +2 o i+3 - coth(2/3) (aia i+1 + aia i+3 )] 

i 

The functions Vk, k — 0, 1, ... contain the interaction terms that arise in fcth order of perturbation 
theory. The system size N enters only through the summation over i. Clearly, the functions Vk 
respect the basic symmetries of our model: the usual Ising Zjzf invariance, and oddness in d under the 
exchange of even and odd sub-lattices. Of course, Vq is the Ising Hamiltonian. In stark contrast to 
the equilibrium case, the stationary distribution q({cr}) cannot be written as exp (— f3H({a})) with 
an appropriate Hamiltonian H({a}), since the potentials Vk possess a complicated dependence on the 
two inverse temperatures (i e and f3 . We also note that, at least for fc < 2, the potentials 14 are sums 
of spin operators each of which involves only a local "cluster" Sk = {o~i, o~i+i, . . . Oi+k+i} of (at most) 
k + 2 spins. Hence, it appears that the size of these clusters grows linearly with k, but the coupling 
strength between two spins, Oi and o~i+k+i, decreases exponentially as d k . We conjecture that this 
behaviour lies at the source of the exponential decay which characterises correlation functions in this 
model gHU. 

We conclude this section with some comments on the accuracy of our perturbation theory. The 
effective expansion parameter, up to and including second order, is Ad, so that the mean inverse 
temperature, (3, also plays a role. In particular, the limit p — > 00 looks troubling at first sight, since 
A oc sinh(4^) (cf. Eqn. (Jlq)). However, a more careful analysis shows that |Ad| < 1 independently of 
how 7 and 7 e approach unity, and our perturbation theory remains valid for low mean temperature. 
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FIG. 1. Quantity \, from Eqn. (|2l|), evaluated for a chain of 6 sites for first (diamonds) and second (circles) 
order perturbation theory. We plot \ as a function of d, at fixed mean temperature 7 = 0.4. 

To test its numerical accuracy, we compare our perturbative solution, g per t{{o'}), with an exact re- 
sult, q ex ({a}), for a short chain of 6 sites | 2^| . For such small systems, the zero eigenvector (stationary 
distribution) of the master equation is easily found. As a quantitative measure, we define 



Keeping the mean temperature, i.e. (j e + 7o)/2, fixed, we expect \ t° scale as \ ~ d k+1 for d — ► 0, 
if we consider q per t. to kth order. Figure |l| confirms this scaling behaviour for k = 1 and k = 2, at 
(7e + 7o)/2 = 0.4, for a large range of d. Remarkably, even for the largest d investigated (0.4), x 
remains well below 0.01, indicating that our perturbative solution is an excellent approximation to 
the exact result. For smaller d, the deviations are even smaller, giving us considerable confidence in 
our perturbative approach. 

In the next two sections, we turn to two general features of our two-temperature spin chain, namely 
the violation of detailed balance and the concept of entropy production. 



Following [25f| . it is straightforward to show that our dynamics violates detailed balance, even 
without knowing the steady-state distribution. The proof involves closed orbits (loops) in configuration 
space, {<t}i — > {o-} 2 —►■••—> {(r}k — * of "length" k = 1,2, ... Each of the configurations in the 

loop differs from its successor by a single spin flip. For each closed orbit we consider the product of 
the transition rates, traversing the loop in the "forward" and "reverse" direction, namely 




(21) 



IV. VIOLATION OF DETAILED BALANCE 
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fc-1 

Uf = w ({a} k -> {ajx) J] ui ({<r}, {a} j+1 ) 

3 = 1 
fc-1 

w (Mi M*0 J| w (Mi+i Mi) 



n,. 



respectively. Detailed balance is violated, if we can identify a single loop for which 11/ 7^ II r . For 
our model, the shortest loop of this type has length 4 and involves four nearest-neighbour spins, 
(Ji—i, ...<Ji+2- We begin with a configuration {cr}i in which these four spins are all +1, and i is even. 
The two central spins play the key role: {<r}2 has spin Ui = — 1; {a}s has at = —1 and er,;+i = — 1; 
finally, {a} a has a i+1 = — 1. Flipping cr i+ i back to +1 restores In defect language, we first 

create an even defect pair, then the right defect, followed by the left one, diffuses right by one lattice 
site, and then the resulting odd defect pair is annihilated: 000 — > 110 — > 101 — > 011 — > 000. It is easy 
to compute 11/ = (1 — 7 e )(l + Jo) which differs from n r = (1 + 7 e )(l — 7o). Thus our model does 
indeed violate detailed balance. 

Our (perturbative) steady state solution sheds further light on this issue, demonstrating which, 
and how, elementary spin flips contribute. We focus on the net probability current F between two 
configurations, {a} and {crW}, 



(22) 



In zeroth order of perturbation theory, all F vanish since detailed balance holds. To obtain the leading 
non-vanishing contributions to F, we decompose the rates according to Eqn. (^J) and insert the first 
order expression for the stationary density. To 0(d), we find 



F 



(Vh^ 1 } 



with 



Z~ 1 d(A 1 +A 2 )q Q ({a})+0(d 2 ) 

?o(M ] }) 



A\ = Aj((7j_l, <Tj, <7j + l) — Aj(0i_l, — CTj, CT, + i) 
A 2 = U)((Ti_l,(Ti,(Ti + l)Ql ({tr}) - U)(o-i_l 



90 (M) ' 

9o ({*M}) 



Cj, 0"j+l) 



?o(M) 



Qi(W li] }) 



Here, Ai reflects the effect of the non-equilibrium rate, A, whereas A2 contains the first-order correc- 
tion, Qi ({<r}), to the Ising stationary distribution. Using our definitions and results from the previous 
section, one finds easily that 



(-1)' 



o~i + 0i+i) [l + exp (-2/3(7, (<7i_i + Oi +1 ))] 



(23) 



and 



A, 



sinh(4/3) 
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(o"i-i + fj+i) 



{ [n ({a®}) - n e ({a®})] - MM) - n e ({a})}} 



(24) 



One can show, by considering elementary spin flips and their effect on the defect representation, that 
the {...} bracket in Eqn. (24) can take only the values ±1 or 0. 

We are now ready to evaluate the relevance of different spin flips types (cf. the discussion following 
Eqn. (||)) for detailed balance violation. Spin flips of type (a) and (b) are inverse to one another: 
where (b) generates a defect pair, (a) annihilates it. To be specific, we let {tr} — > {<r^} be a type (b) 
transition, with = <Ji = o"i+i in {a}. For this configuration, we find 



A 1 = (-l) i [l + exp (-4,9)] , 

A 2 = \ { [n ({c-®}) - n e ({c-W})] - MM) " «e(M)]} 



(25) 
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Given the constraint on the {...} bracket, the inequality | ^4_2 1 < 2 |^4i| follows. Thus, A\ determines 
the sign of the probability current. For T e > T a (d > 0) we have F ({a}, {ct' 1 '}) > if i is even, and 
negative otherwise: At the hotter (even) sites, defect pairs are more often generated than removed. 
At the cooler (odd) sites, the situation is reversed. In this manner, spin flips of types (a) and (b) 
always contribute to detailed balance violation in the presence of two temperatures. 

If the spin flip {a} — > {c' 1 '} is a type (c) process (i.e., defect diffusion ), then its inverse also belongs 
to this class, and both occur with rate 1. Since cx^-i = — for a type (c) spin flip, A\ vanishes 
and only Ai contributes to the probability current. As a result, a type (c) transition violates detailed 
balance in first order of d (i.e., generates a nonzero F) only if it changes the difference of odd and even 
defect pairs. In defect notation, the process ...0, 1, 1, 0, 0, ... — > ...0, 1, 0, 1, 0, ...violates detailed balance 
while ...0, 1, 0, 0, ... — * ...0, 0, 1, 0, ... does not. If the defect pair in the former process is centred on an 
even site, we will have F ({a}, {cr^}) > for T e > T a , etc. 

Including terms of 0{d 2 ), we find that all type (c) spin flips now generate non- vanishing probability 
currents. Thus, currents between pairs of configurations, connected by single spin flips, are generically 
nonzero. 



V. ENTROPY PRODUCTION 

We finally consider two different approaches to entropy production in our model. The first p[ is 
very intuitive and tailored to our specific model, by focusing on the energy ( "heat" ) flux from one 
bath to the other. The second method |^,^ M is more general, defining an entropy production for 
any master equation. In the following, we show that both approaches are equivalent for our model. 

We briefly review the first method. As the spin chain is connected to two baths at different temper- 
atures, heat flows from the warmer reservoir into the spin chain and on into the colder one. Since the 
energy difference of two configurations connected via a single spin flip at site i is —2 Oi [pi—\ + <Ji+i), 
the energy transfer at this site is given by 

= (-2<7i ((Tj-i + <T l+ i) w i (a l -> -a;) ) (26) 

where (•) denotes the configurational average with respect to q({o~}). Standard thermodynamics then 
suggests to define the entropy production of the whole chain as 

. = N(n^ + i) + m\ . (27) 



2 V T T e 

Using Eqn. (Q) for the rates, is easily expressed in terms of the exactly known [ fl4[ two spin 
correlation function 



(<j i <r j ) = y/A(i) Ajj) 

with A(4\ = / ^ = (7e+7o)/(2 7o), i even 
with A{%)- | Ao = (7e+7o)/(27e)i i odd , 

A = -L= (l - yfl - le lo) • (28) 



Hence, 



*-t<* -*>(£-£) (29) 

follows exactly, for any choice of temperatures. We note that the entropy production is positive 
whenever the two temperatures differ. 

The more general approach [|l6];^6] |||] to entropy production begins with the Gibbs entropy for the 
full time-dependent probability distribution p ({cr}, t): 

S G {t) = -J2p(Wht)lap({a},t) . (30) 



For the stationary distribution q ({cr}), the time derivative dSc/dt must vanish. With the help of the 
master equation, we find dSc /dt = Si — S% where 



Si 



x In 



Wi(-ai -* Oi) g({crW}) 



The two expressions Si and S2 constitute two alternate representations for the entropy production. 

In particular, we recognise the probability current F({a}, {c^}) of Eq. ( ^ ) in the expression for Si- 

One easily sees Si > since each term in the sum is non- negative. Moreover, Si vanishes if and only 
if detailed balance is satisfied. Thus, strictly positive entropy production is equivalent to the violation 
of detailed balance. 

We now evaluate the expression for $2 in Eqn. ( |3l| ) p9| . Since our model exhibits translation 
invariance with period 2, we only have to consider one even (2z) and one odd (2i + 1) site, so that 

^2 = y^ Yl qda^Wjiaj -> -aj) [-2/3 J -£T j (cr i _ 1 + a j+1 )} . (32) 

{cr} je{2i,2i+l} 

Here, we have used the simple relation 

Wi{<7i -> -CTj) 

— = exp -2fta-i (cj_i + cr i+ i) 

Wi( — (Ti — > CTj) 

Recalling Eqn. (|2^), we see immediately that 

- iV 

S2 = y [/3 e i(2i) + /3 i(2i+l)] (33) 

which is exactly the expression (ETJ) that was derived from thermodynamics. Thus, we confirm the 
consistency of the more abstract general expressions ([H]) and the intuitive thermodynamic approach. 
While this enhances our confidence in (|3l|), a better understanding in terms of more basic physical 
principles is still outstanding. 



VI. CONCLUSIONS 

To summarise, we have gained some analytic insight into the steady state configurational proba- 
bilities of an Ising spin chain driven out of equilibrium by a coupling to two heat baths: Subject to 
different temperatures, spins on even and odd sites are updated according to a generalisation of the 
usual Glauber rates. A perturbative calculation in the temperature difference shows that the station- 
ary distribution for this non-equilibrium model is rather complex, with longer-range and higher-order 
spin operators appearing. Clearly, the energy of a configuration no longer determines its statistical 
weight. Unfortunately, at this stage we can only conjecture the structure of this distribution beyond 
second order in perturbation theory. It seems likely that, at each order, additional spin operators 
will have to be introduced. Turning to other characteristics of this non-equilibrium steady state, we 
could understand in more detail how the presence of two temperatures violates detailed balance. We 
could also show the equivalence of two unrelated definitions of entropy production, thus providing 
additional support for the general expression Jl6| , p6| -p8| . 

Of course, numerous questions remain open. In particular, we would like to reconcile the apparent 
proliferation of longer-range spin couplings in the configurational probabilities with the observed 



"triviality" of long- wavelength properties. At present, we have strong indications |23| that arbitrary 
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spin correlation functions decay exponentially, quite similar to the equilibrium Ising chain. However, 
the most fundamental question still remains unanswered: are there any a priori criteria, beyond simple 
symmetries, that determine which configurations will appear with equal weights? In equilibrium, the 
determining quantity is their energy. Far from equilibrium, the key concept is still missing. 
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